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Abstract. In this contribution, we discuss the finite-element discrete variable representa- 
tion (FE-DVR) of the nonequilibrium Green's function and its implications on the descrip- 
tion of strongly inhomogeneous quantum systems. In detail, we show that the complemen- 
tary features of FEs and the DVR allows for a notably more efficient solution of the two-time 
Schwinger/Keldysh/Kadanoff-Baym equations compared to a general basis approach. Particu- 
larly, the use of the FE-DVR leads to an essential speedup in computing the self-energies. 

As atomic and molecular examples we consider the He atom and the linear version of 
in one spatial dimension. For these closed-shell models we, in Hartree-Fock and second 
Born approximation, compute the ground-state properties and compare with the exact findings 
obtained from the solution of the few-particle time-dependent Schrodinger equation. 



1. Introduction 

In the last decade, the application of the nonequilibrium Green's function (NEGF) to describe 
strongly inhomogeneous quantum systems has started to become an actively considered subject. 
Thereby, various finite and localized systems have challenged attention and different many- 
body approximations have been applied. Recent state-of-the-art approaches discuss small atoms 
and molecules [U |2] , few-electron quantum dots [H 0] and quantum dots coupled to leads [5] , 
molecular junctions [6], and Hubbard nanoclusters [7|. In part, these works also include the 
monitoring of the system's temporal evolution which, accounting for correlation and memory 
effects, requires to solve the two-time Schwinger/Keldysh/Kadanoff-Baym equations [8j EJ ITU] 
(SKKBE). 

To solve the SKKBE for homogeneous quantum systems p^l fT2l fT3l fl^l fT5l fT6l fTTl fTBl 
\19\ |20| |2T| [22] has become routine. However, this still does not hold for finite, localized and 
inhomogeneous systems. The reason is, that in contrast to homogeneous systems where one 
spatial coordinate or momentum drops out of the NEGF, the inhomogeneity claims adequate 
resolution in both coordinates or momenta. To meet these requirements all above mentioned 
works rely on an expansion of the NEGF in one-particle orbitals: for atomic and molecular 
systems, linear combinations of Slater-type or Gauss-type orbitals are being used (also in 
connection with tight-binding models), whereas for other classes of problems e.g. potential- 
eigenstate basis functions are being utilized. Beyond a basis ansatz, alternatives such as 
direct grid (finite-difference) methods do not exist as they are computationally very expensive. 



Nevertheless, also basis approaches are so far limited to relatively small basis sets. Generally, 
the numerical complexity involved in the description of the binary interactions does not permit 
to extend nonequilibrium calculations to much larger basis dimensions than a guide number of 
rib < 50 orbitals. This explains the low resolution in the description of photoionization processes 
of model atoms |23| YM\ using nonequilibrium Green's functions. 

In this contribution, we — for the first time — develop a grid-based approach in the frame of 
the finite-element discrete variable representation (FE-DVR), e.g. Refs. [25j [261 113 I2H] and 
references therein. This method leads to specially-designed, flexible basis sets which are capable 
to combine the advantages of pure grid and standard basis approaches, see Sec. [2] and also 
Ref. [29]. In particular, it allows for a very efficient treatment of the binary interactions and, in 
turn, a drastic simplification of higher-order self-energy expressions, which, generally, require the 
main effort in all NEGF calculations. As a result, only 0(n%) semi-analytical matrix elements of 
the interaction energy operator are required to compose the second Born self-energy. This has to 
be compared to a general basis representation of the NEGF: There, O(n^) matrix elements are 
involved, cf. [lj. As a consequence, the use of the FE-DVR enables more efficient calculations 
at less storage memory and computing time and provides the basis to consider also spatially 
extended hamiltonians, where particles may occupy large domains in coordinate space. 

After the outline of the method, we, in Sec. [3j apply the FE-DVR representation of the NEGF 
to compute the ground-state properties of atomic and molecular models. First, in Sec. 3.1, we 



discuss the one-dimensional (ID) helium atom [301 ED E2 E31 EH E3 ES EZ] , where we focus on 
technical details such as grid size, DVR basis size, and convergence in the case of the Hartree- 
Fock and second Born approximation. The benchmarking results shown are also of relevance 
for time-dependent calculations, as they define and border the requirements to resolve explicit 
correlation effects (within a NEGF approach) e.g. the two-electron resonances [38] which are 
embedded in the one-electron continuum of the atom and are going to be occupied during laser- 
atom interactions — see also [39] in the present volume. Finally, we consider the linear molecular 
ion in the symmetric ID singlet configuration [401 W\\ as an example of a two-electron 
molecule and vary the interatomic distance to record the binding-energy curves in Hartree-Fock 



and second Born approximation, see Sec. 3.2 From this we can extract the minimum ground- 
state energies and the respective bond- lengths, which are compared to the findings from the 
few-particle time-dependent Schrodinger equation. 



2. The nonequilibrium Green's function in FE-DVR representation 

For the description of the NEGF using finite elements together with the discrete variable 
representation [32] (DVR), we consider the general iV-electron Hamiltonian 

N 

i=l i<j 

with the kinetic energy ti = — Vf/2, the time-dependent potential energy V{ = v(xi,t), and 
the pair interaction energy Uij = u{\xi — Xj\) formulated in atomic units. The one-particle 
nonequilibrium Green's function with space-time arguments 1 = (x,t) and 1' = (x',t') reads 

G(l,l') = -i (^(1)^(1')) , (2) 

where spin is omitted, ip(l) and are electron field operators, and Tq denotes time-ordering 

on the full Keldysh contour C. According to system ([l), G(l, I') obeys the SKKBE [5] El [TO] 



{id t -H(l)} G(1,1') = 5 C (1-1') + jfd2E[G](l,2)G(2,l') , 



(3) 



X" 
• 



= n e n q - 1 



at 



Figure 1. Interval [0, xq], as discretized in FE-DVR representation with a number of n e finite 
elements and n g generalized Gauss-Lobatto points (per element), denotes the dimensionality 
of the extended basis covering the whole interval. 



with addition of its adjoint equation. Further, H{\) is the one-electron (kinetic plus potential) 
energy, E(l, 1') denotes the self-energy, and equilibrium initial correlations are treat in the mixed 
Green's function approach [131 SH US] , cf. Sec. |2.3| 

The favorable aspects of the FE-DVR representation have been successfully used for the 
time-dependent Schrodinger equation (TDSE) by Rescigno [25] and others, e.g. [28]. There, the 
accuracy of the DVR [46 on the one hand and the sparse character of FEs on the other led to 
an efficient TDSE code which is well parallelizable. In our case, this hybrid approach allows us 
to rewr ite the SKKBE @ in a highly effective matrix notation using optimal combinations of a 
grid and a local basis. 



2.1. Grid-based ansatz 

On a predefined spatial ID interval I = [0,xo], we expand the nonequilibrium Green's function 
of hamiltonian ([!]) as 

G 0» = E E AO*) xZM) > x, x 'ei, (4) 

i 1 m 1 i 2 m2 

with time-dependent complex coefficients <7™ m /(i, t') and real basis functions Xm( x )- Outside the 
interval I, we assume the NEGF to vanish. All indices i (superscripts), ranging 0, 1, . . . , n e — 1 
in Eq. Q, are linked to a grid of length xo composed of n e finite elements, see Fig. [T] All 
indices m (subscripts), ranging 0, 1, . . . , n g — 1, are connected to locally defined basis functions 
Xm( x ) wn ich are constructed as follows: First, we divide the interval X into finite elements with 
boundaries x° = < x 1 < x 2 < . . . < x Ue ~ 1 ,x ne = xq. To each FE [x*, x l+1 ], we then attach a 
local DVR basis using the generalized Gauss-Lobatto (GGL) points [25] x l m and weights w l m : 

l -{{x i+l -x i )x m +{x i+l +x 1 )} , (5) 

f - x*) , (6) 

with the standard Gauss-Lobatto points x m (and weights w m ). For the special case of Legendre 
interpolating functions, the points x m are defined as roots of the first derivative of Legendre 
polynomials P n {x) according to 

= , (7) 



and the weights w m are given by 



where n g denotes the total number of basis functions per element. In our approach, we keep the 
number of basis functions constant in each FE |47] ■ However, in order to combine the locally 
defined DVR basis functions to a continuous basis set, the FE-DVR space is being spanned by 
two classes of functions — 'bridge' and 'element' functions, see also Ref. [29 . The bridge function 
(the case m = in Eq. Q) extends over two adjacent FEs and ensures communication between 
the grid domains i and i + 1. In particular, it guarantees the continuity of the NEGF. The 
'element' functions are zero at and outside the respective element boundaries. Using definitions 
(|5j and (|6]), the basis functions have the explicit form 
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and are orthonormal in the sense of the generalized Gauss-Lobatto quadrature. In Eq. ([9]), the 
so-called Lobatto shape functions [23 HE] fm( x ) are defined as 



fm( x ) 
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x < x l and x > x 
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and obey fmi x m') = feW' As in the last finite element n e — l (i.e. in the element [x n f~ 1 , x ]) 
no bridge function is needed due to the boundary condition of vanishing (7(1, 1') outside the 
interval /, the total FE-DVR set consists of 
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basis functions. Finally, we note, that a generalization of Eq. ([9]) to higher dimensions (2D and 
3D) is, in principle, possible by using a product ansatz for the coordinate functions [27]. 

2.2. Matrix elements of the kinetic, potential and interaction energy operators 
With representation Q of the NEGF, the SKKBE will tran sform into an equation of motion 

Obviously, this equation involves 



2.3 



for the matrix g(M') with elements gm? m2 (t,t?), cf. Sec. 
also the kinetic, potential and interaction energy operator of Eq. ([T]) in matrix form, which we 
specify in the following. Thereby, integrations over coordinate space are performed by using 
the generalized Gauss-Lobatto quadrature rule, and case differentiations arise from the basis 
functions Xm( x ) being split into element and bridge functions. 

First, let us consider the potential and the kinetic energy in FE-DVR representation: The 
potential-energy matrix is given by 
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with 
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This implies that the potential energy is diagonal with respect to elements i and local DVR basis 
indices m, and that, consequently, it can be represented by a vector of dimension rn,. Moreover, 
Eq. (12) holds true also for any other local operator. As the operator of the kinetic energy 
is non-local in coordinate space, the matrix elements t l ^ m are not diagonal. We follow the 
derivation of Ref. [25] and obtain the kinetic-energy matrix as 



fW2 



■rii 



- /)• • f* 1 111!* 1 7I1* 1 

2 u *i»2 b ra\m2 L mi m 



-1/2 



"12 J 



2 I °*l«2 fc n„-l,m2 ~ r U U«2-1 L 0m, 2 



W 



tin — 1 



+ w 



-1/2 



2 ( ^1*2 dn 3 -l + ^1*2 + 1 Ci( 



mi 



<-l + 



-1/2 



8 ri 



'-' I *ri5-l,n 9 -l + * 



00 



+ ^ui 2 -l *oJn g -l + ^U«2 + l *ng-l,0 



77', 



+ W. 



h+1 



W, 



+ w% +1 



1/2 



(14) 

m\ > 0, 7722 > 
772-1 = 0, 777,2 > 
777i > 0, 7772 = 

, 777i = "2,2 = 



Here, the case differentiations lead to the matrix having a block-diagonal form [27., and the 
quantity i % mim2 is given by 



"mim2 



dx 



dx 



(15) 



which involves the first derivatives of the Lobatto shape functions at the GGL points, see also 
Ref. [25]. 

Next, let us focus on the matrix elements of the binary-interaction operator u (the two- 
electron integrals) which are carrying a set of four index-pairs (i, m) and are defined by 



Hl2,l3U 
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(16) 



In a general basis approach, the two-electron integrals [Eq. (16) with all index-pairs replaced 
by single indices] often require a careful analysis, as they are not analytically accessible and 
have to be numerically precomputed for all combination of indices, e.g. [3j. Although symmetry 
relations (39] help to restrict oneself to a subset of indices, the effort scales with 0(n^) and 
thus can be huge for larger basis sets. On the contrary, using the FE-DVR, the evaluation of 
the two-electron integrals turns out to be much simpler. In particular, the integrals can be 
performed in a semi-analytical way such that Eq. (16) reduces to 

,n«2,«3*4 



u 



mi rri2, 0237714 



m,2 u m3m4,U mim2 ■ 



«1*2 



where the kernel matrix u is symmetric and follows as 
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To obtain Eq. (18), we have used the separable form of the discretized interaction potential 
u(\x — x'\), and, correspondingly, the quantities a l m denote the eigenvalues of the matrix 
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and i are related to the eigenvectors f3 % 



6' 



ai i 
• vn.'i 



, m > 



hi i 



m'{n g ~l) W n g 



-1 + Pm'O 



i '(*+ 1 )„„«+l 







w: 



.! + W\ 



i+l 



(20) 



m 



In comparison with any single-electron matrix element (such as the kinetic or potential energy) , 
in FE-DVR, the calculation of the binary-interaction matrix elements involves just an additional 
but numerically elementary matrix diagonalization. Furthermore, besides the fact that with 



Eq. (17) the two-electron integrals attain a very simple form independent of the specific pair- 



interaction potential, only a single matrix of dimension n\, x needs to be stored in the 
code. This memory-friendly property is based on the high degree of diagonality determined 
by the product of Kronecker deltas and represents a main attractive feature of the FE-DVR 
representation. Particularly, this aspect opens the way towards efficient NEGF calculations, 



since Eq. (17) has direct consequences for the structure of the self-energies, see the following 
Section. 



2.3. Equations of motion 

Once all relevant matrix elements are known with respect to the chosen FE-DVR basis (|9]), 
we can start to solve the Schwinger/Keldysh/Kadanoff-Baym equations Q for the one-particle 
Green's function G(l, 1') expanded in the form of Eq. Q. This implies, though, the SKKBE in 
the matrix form of the finite-element discrete variable representation: 
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where we have denoted 6*^ , 

with its adjoint equation. Further, Eq. (22) separates the self-energy matrix £^ 2 m2 (£, t') into a 



h(t) = t + v(i), and Eq. (21) has to be supplied 



time-local Hartree-Fock part (X! HF ) and a contribution 5] corr that accounts for electron-electron 
(e-e) correlation and memory effects. However, as an exact treatment of e-e correlations is 
impractical, we have to apply many-body approximations for which the second Born diagrams 
provides one of the most basic models. Hence, besides the general form of the HF self-energy 
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where t + denotes t — > t + e>o we, in second Born approximation, have 
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Eqs. (|23| and (24) involve the spin degeneracy factor a £ {1,2}, the matrix elements u l ^ 2 m2 of 
Eq. (|17p and show a very simple form compared to the situation when a general basis is applied, 
e.g. jl]. The reason for this is the subtle structure of the FE-DVR. In detail, the Hartree 
term is completely diagonal [including a single sum over elements] and the exchange term 



involves only a product of two matrix elements. For the second-order Born terms, the degree of 
simplification is even more drastic: In a general basis representation, two sums are required for 
each full vertex point in the second-order self-energy diagrams and, additionally, a single sum 
is needed for the start- and the end-point. This leads to a scaling with 0{nT). However, in our 



case, due to the diagonality of the two-electron integrals, cf. Eq. (17), the evaluation reduces 
remarkably to a scaling with 0(n 2 ) per matrix element. 

For the atomic and molecular model calculations to be outlined as first benchmarks of the 
method in Sec. [5J we, in this contribution, restrict ourselves to the (equilibrium) ground-state 
properties. For completeness, we give a short summary of the main computational steps. First, 
in the FE-DVR picture, the Hartree-Fock equilibrium Green's function, denoted g°(r), follows 
from 

^m^mafe (0 )] = ^mi"i2 v mim 2 ^^mimafe (0 )] i (25) 

where gmim 2 { T ) is the HF approximation of g^' 2 mo (t, t') \ t _ t i =0+iT with r G [— /?, 0] and (3 is the 



inverse temperature, and S° is defined via Eq. (23), but with g being replaced by matrix g°. 



Thereby, Eq. (25) has to be solved to self-consistency by iteration and as result we obtain 
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Here, the vector contains the energy eigenvalues of h , the matrix c^ m , summarizes the 
corresponding eigenvectors, and the chemical potential \i is determined by normalization of the 



Fermi distribution: N = ^2 im j "$(An — A 4 )- Eq. (26) solves a simple differential equation, see 
e.g. [3], and corrections due to e-e correlations in second Born approximation are obtained by 
insertion into the Dyson equation [33] (the SKKBE in the limit t — t' = + ir) for the full 
Matsubara Green's function g M (r): 

M,hi 2 ( \ _ 0,iii2 ( ) , T(2)hh( ) (07) 

{■13 

j{1)*m( t \ = / dfa°' lul (T-f)I Wll2 (f) 



im 







lWiih( T ) _ / df Iy, m ' 1i1 (t - f) - S(t - f) S 0,Jl * \ a M ' U2 (f) 



where S M (r) = 5(r)I] HF (r) + E corr (r) with g M instead of g in Eqs. g and @. The 
convolution integrals in Eq, (27) are performed in sequence by direct integration. This has been 
found to be more stable and controllable than the method applied before in Refs. [21 Hj. Once 
the self-consistent g M (r) is computed, we have direct access to many observables, e.g. the one- 
electron density is obtained as n(x) = ^2 hmi E; 2 m 2 X mi ( x ) Xm,2 ( x )9m{m2 (0 )• The total energy 
is computed similar as in Refs. |43| H] . Overall, in order to ensure the atoms and molecules being 
in the ground state, we set the inverse temperature (5 = 100. 



3. Performance for atomic and molecular models 

As first benchmarks and preparatory work for the investigation of the temporal evolution of 
small atoms and/or molecules following an external (e.g. laser- induced) perturbation, we here 
consider their equilibrium-(initial-)state preparation within the FE-DVR context. As examples, 
we focus on the He atom and the linear molecular ion H^, modeled in one spatial dimension. 
For both two-electron systems, the Coulomb potential is considered in the regularized form 
u(\x — x'\) = \{x — x') 2 + l] -1 / 2 which, from the physical point of view, allows for a transverse 
extension of the few-particle wave function. Furthermore, e-e correlations are treated in second 
Born approximation. 




Figure 2. a) Convergence of the Hartree-Fock ground-state energy for the ID helium atom 
against basis dimension n.&. As indicated, each curve corresponds to a different arrangement of 
the FEs within a total interval of 200 a.u. length. For the upper set of lines the FEs are equally 
distributed, whereas for the lower set the central FE is one atomic unit wide, and the width of 
neighboring elements is linearly increasing towards the boundaries. In addition, the number of 
local DVR basis functions has been varied between n g = 2 and 15, see numbers on the curves, 
b) He ground-state energy in second Born approximation (n^ = 202) with respect to different 
T-grid parameters u and p. 



3.1. The ID helium atom 

The helium atom is the most elementary closed-shell two-electron system. In one spatial 
dimension, it is well modeled by the nucleus potential v(x) = —Z [{x — xq/2) 2 + /j] -1 / 2 , where 
the atomic number is Z = 2, x £ [0, xq], and p is a regularization parameter. In this setup, the ID 
helium atom serves as a fundamental 'testing ground' for multi-electron calculations |36[ [37| [39] 
and provides many features of the single- and double-ionization dynamics including the so- 
called knee structure [34, 35 . Considering the singlet state, we refine the model and set a = 2 
in Eqs. (23) and (24), and used p = 1. 

Further, for the equilibrium calculations, we have used a FE-DVR basis that covers a total 
interval of 200 a.u. length. This corresponds to a domain that is about 100 times larger than 
the characteristic extension of the ground-state wave function or density, cf. Ref. Such 
a grid extension is more than adequate to resolve the ground-state features of the model (to 
be discussed here) but will become crucial, when the helium atom is perturbed by external 
fields and, in turn, electrons start to occupy highly excited or continuum states. In this sense, 
our results are benchmarks also with relevance for the computation of the system's temporal 
evolution. In particular, a grid with an extension of about 200 a.u. should be well capable to 



resolve the two-electron resonance states embedded within the one-electron continuum of dipole 
spectra. This follows from the performance of the few-particle time-dependent Schrodinger 
equation (TDSE) using absorbing potentials, e.g. [50 and references therein. We note, that also 
in the FE-DVR approach, such imaginary one-electron potentials that damp reflections at the 



interval boundaries are easily implemented, just allowing the matrix elements in Eq. (12 ) 

to be complex. 

The explicit partitioning of the interval into finite elements has been organized as follows: 
(case I) the interval is divided into equidistant segments, (case II) the central FE has a width 
of one atomic unit and the width of the surrounding elements is linearly increasing towards the 
interval boundaries. The effect of these segmentations on the Hartree-Fock ground-state energy 
convergence of the helium atom is displayed in Fig. [2] a) as function of the total basis size. In 
case I, for equidistant FEs, the convergence is relatively slow with and strongly depends on 
the number of elements as well as on the number of local DVR basis functions used (see the 
different symbols and lines). Consequently, more than > 550 functions are needed for the HF 
energy to deviate by less than 10 -4 Ha from the converg ed HF result JSg£ t = -2.2242096 Ha. 
In case II, the situation is completely different as essentially more basis functions are available 
in the center region of the interval. This enables a superior representation of the Matsubara 
Green's function G m (x,x';t) in coordinate space, and leads to adequate convergence at 200 to 
300 FE-DVR basis functions with an error reduced by several orders of magnitude compared 
to case I. Furthermore, the ground-state energy depends less on the number of elements. For 
completeness, we note that, besides the total energy, an additional indicator for the basis quality, 
is to look at how well the potential v (x) can be expanded into the chosen FE-DVR basis. 

For case II with n e = 29 and n g = 7, we have computed the ground-state energy in second 
Born approximation. Thereby, all quantities are found to be well converged with respect to the 
basis size, compare with the HF case in Fig.[2]a). The time -argument r in the Matsubara Green's 
function (r S [— (3,0]) has been discretized using a uniform power mesh [51] with parameters u 
and p — for definition see e.g. Refs. [U 143] . Fig. [2]b) indicates the convergence with respect to 
these parameters, where the total number of r-grid points is given by 2up +1. At j3 = 100, a 
mesh parameter p > 10 ensures the particle number N = ^2 im grnm(0~) being sufficiently stable 



during iteration of the Dyson equation (27). Particularly, with more than 1000 grid points it is 
possible to compute the ground-state energy to relatively high precision, E 2ndB = —2.233419 Ha. 
As result, the electron-electron correlations lower the total energy accounting for 66% of the 
correlation energy and, hence (improving the HF result), approaches the exact ground-state 
energy which is —2.2382578 Ha. For the discussion of other observables such as the one-electron 
density in HF and second Born approximation, the reader is referred to Ref. [29 . 

3.2. The linear molecule H% 

With more than two nuclei, the molecular ion H3" can, in ID, only be realized in its linear 
version [52] , where the H-H bonds are oriented parallel to each other. Hence, the electrons move 
along the molecular axis and the one-electron potential is modeled as 
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v d (x) = [(x - (x - d)/2) 2 + lj + [(x- (x Q + d)/2) 2 + lj + [x 2 + 1]" 1/2 + - , (28) 

where d/2 denotes the interatomic distance, and the last term (offset) collects all internuclei 
interaction energy contributions. For the ground-state NEGF calculations, the coordinate space 
has been constrained to an interval of 50 a.u. with a grid of n e = 13 FE being linearly increasing, 
starting from a one atomic unit wide central element. In total, rib = 142 FE-DVR basis functions 
have been used. 

The total binding energy E\, of the singlet state is shown in Fig. [3] a) against distance d 
for the exact solution of the TDSE (dotted [and triangles]), the Hartree-Fock (dashed) and the 




Figure 3. a) Binding-energy curves for the linear Hg~ molecule in one spatial dimension for 
different approximations. The dotted line corresponds to the exact solution obtained from the 
few-particle time-dependent Schrodinger equation (TDSE). The triangles denote the binding 
energy of Ref. [ID], b) One-electron density at the equilibrium bond-lengths 2d\>: 4.3654 (HF), 
4.5579 (second Born), and 4.7698 (exact). The dots indicate the self-consistent positions of the 
nuclei. The gray curves indicate the corresponding potential Vd(x) [scaled by 0.35 and shifted]. 



second Born approximation (solid). Over a broad range of internuclear distances, the second 
Born approximation, thereby, accounts for about 60-70% of the correlation energy, and — leading 
overall to a larger bond — essentially improves the HF result. 

Furthermore, with a value of —1.3396 Ha, Hg~ has the same dissociation threshold as the 
ID hydrogen molecule [29], but, in addition, leaves behind a positively charged hydrogen ion. 
However, we note, that the Hartree-Fock and the second Born approximation, cannot resolve 
this threshold. The reason is that the H^j~ molecule dissociates into open-shell fragments — two 
hydrogen atoms and a single hydrogen ion. These cannot be represented within a spin-restricted 
calculation (with a = 2), and, hence, lead to a strong deviation of E^ in the limit of large d. 
Regardless of this failure of the ansatz, the NEGF calculations are well capable to describe the 
behavior of E\> around the minima in the binding energy curves. Also, the equilibrium positions 
are consistently reproduced, and the correct trend is observed when e-e correlations are being 
included: the bond-length d\y shifts to larger nuclear separations, for the specific values obtained 
see caption of Fig. [3] 

The electron ground-state density of the molecular ion is displayed in Fig. |3]b) together with 
schematic curves for the spatial potentials Vd(x), where d is twice the equilibrium internuclear 
distance, cf. Eq. (28). In HF approximation (dashed curve), the density shows a pronounced 
maximum in the central region of the molecule, whereas the exact density (dotted curve) is 
essentially less peaked. However, also the exact result does not indicate onset of electron 
localization, i.e. does not show separated maxima in n H +(x). In second Born approximation, 
corresponding to a lower total energy (cf. Fig. [3] a)), we resolve the correct trend of this density 



reduction. Moreover, replacing the self-consistent bond-lengths by the corresponding length 
obtained from the TDSE (dash-dotted curves) only slightly improves the results for the Hartree- 
Fock and the second Born approximation. This explains that the substantial differences in the 
density profiles are unambiguously due to e-e correlation effects. 



4. Conclusion 

The FE-DVR ansatz Q provid es an elegant and very efficient way to treat binary interactions in 
NEGF calculations for inhomogeneous quantum systems. To this end, the method uses a flexible 
combination of grid (FE) and basis (DVR) strategies, which allows for simple, (semi)-analytical 
formulas for the matrix elements of the kinetic-, potential- and, especially, the interaction-energy 



operator, cf. Sec. 2.2 Further, for the most basic model that accounts for particle-particle 
correlations — the second Born approximation — the use of a FE-DVR basis enables remarkable 
scaling properties: Only 0{nfy summations are required for the computation of a single matrix 
element of the second-order self-energy instead of 0(n®) summations that are needed in a general 
basis approach. 

Also, we emphasize that the FE-DVR space can, e.g. via finite-element variations, be well 
adjusted to the problem considered and, thus, allows for an efficient expansion of the NEGF 
[also for spatially extended hamiltonians] and quick convergence of the observable s of inte rest. 
This has been exemplified for the He atom and the triatomic molecule in Sec 



3.1 



and 



3.2 



In addition, the method is found to be stable also for large grids [large basis sets with raj, > 500] 
and also for larger particle numbers N < 20, considering interacting fermions in a harmonic 
trap potential. 

Finally, we believe that the FE-DVR method is attractive also for other classes of many- 
body approximations, such as GW or T-matrix calculations, as it will likewise simplify the 
computation of self-energy contributions of higher than second order. Moreover, though we, in 
this contribution, focused on the solution of the Dyson equation, the formalism presented is, in 
particular, well applicable in nonequilibrium situations solving the full two-time SKKBEs which 
will be demonstrated in a forthcoming publication. 



References 

[1] A. Stan etal, Europhys. Lett. 76, 298 (2006). 

[2] N.E. Dahlen, and R. van Leeuwen, Phys. Rev. Lett. 98, 153004 (2007). 

[3] K. Balzer, and M. Bonitz, J. Phys. A: Math. Theor. 42, 214020 (2009). 

[4] K. Balzer, M. Bonitz, R. van Leeuwen, A. Stan, and N.E. Dahlen, Phys. Rev. B 79, 245306 (2009). 

[5] P. Myohanen, A. Stan, G. Stefanucci, and R. van Leeuwen, Europhys. Lett. 84, 67001 (2008). 

[6] K.S. Thygesen, Phys. Rev. Lett. 100, 166804 (2008). 

[7] M. Puig von Friesen, C. Verdozzi, and C.-O. Almbladh, Phys. Rev. Lett. 103, 245306 (2009). 

[8] PC. Martin, and J. Schwinger, Phys. Rev. 115, 1342 (1959). 

[9] L.V. Keldysh, Zh. Eksp. Teor. Fiz. 47, 1515 (1964) [Sov. Phys. JETP 20, 235 (1965)]. 

[10] L.P. Kadanoff, and G. Baym, Quantum Statistical Mechanics (Benjamin, Inc., New York, 1962). 

[11] P. Danielewicz, Ann. Phys. (N.Y.) 152, 305 (1984). 

[12] H.S. Kohler, Phys. Rev. C 51, 3232 (1995). 

[13] P. Bozek, Phys. Rev. C 56, 1452 (1997). 

[14] N.H. Kwong, and M. Bonitz, Phys. Rev. Lett. 84, 1768 (2000). 

[15] M. Bonitz, D. Kremp, D.C. Scott, R. Binder, W.D. Kraeft, and H.S. Kohler, J. Phys.: Cond. Mat. 8, 6057 
(1996). 

[16] D. Semkat, D. Kremp, and M. Bonitz, Phys. Rev. E 59 1557 (1999). 

[17] D. Semkat, D. Kremp, and M. Bonitz, J. Math, Phys. 41, 7458 (2000). 

[18] R. Binder, H.S. Kohler, and M. Bonitz, Phys. Rev. B 55, 5110 (1997). 

[19] N.H. Kwong, M. Bonitz, R. Binder, and S. Kohler, phys. stat. sol. (b) 206, 197 (1998). 

[20] W. Schafer, J. Opt. Soc. Am. B 13, 1291 (1996). 

[21] L. Banyai, H. Haug, and P. Gartner, Eur. Phys. J. B 1, 209 (1998). 

[22] P. Gartner, J. Seebeck, and F. Jahnke, Phys. Rev. B 73, 115307 (2006). 

[23] D. Hochstuhl, K. Balzer, S. Bauch, and M. Bonitz, [arXiv:0902.0768"l Physica E, in press (2009). 



[24] M. Bonitz, D. Hochstuhl, S. Bauch, and K. Balzer, arXiv:0909.1964, submitted to Contributions to Plasma 
Physics (2009). 

[25] T.N. Rescigno, and C.W. McCurdy, Phys. Rev. A 62, 032706 (2000). 

[26] L.A. Collins, S. Mazevet, J.D. Kress, B.I. Schneider, and D.L. Feder, Physica Scripta T110, 408 (2004). 
[27] B.I. Schneider, L.A. Collins, and S.X. Hu, Phys. Rev. E 73, 036708 (2006). 

[28] J. Feist, R. Pazourek, S. Nagele, E. Persson, B.I. Schneider, L.A. Collins, and J. Burgdorfer, J. Phys. B: 

At. Mol. Opt. Phys. 42, 134014 ( 2009). 

[29] K. Balzer, S. Bauch, and M. Bonitz, |arXiv:0910.5458l submitted to Phys. Rev. A (2009). 
[30] M.S. Pindzola, D.C. Griffin, and C. Bottcher, Phys. Rev. Lett. 66, 2305 (1991). 

[31] R. Grobe and J.H. Eberly, Phys. Rev. A 48, 4664 (1993); S.L. Haan, R. Grobe, and J.H. Eberly, Phys. Rev. 

A 50, 378 (1994). 
[32] D. Bauer, Phys. Rev. A 56, 3028 (1997). 

[33] W.-C. Liu, J.H. Eberly, S.L. Haan, and R. Grobe, Phys. Rev. Lett. 83, 520 (1999). 
[34] N.E. Dahlen, and R. van Leeuwen, Phys. Rev. A 64, 023405 (2001). 
[35] M. Lein, E.K.U. Gross, and V. Engel, Phys. Rev. Lett. 85, 4707 (2000). 

[36] J. Zanghellini, M. Kitzler, T. Brabec, and A. Scrinzi, J. Phys. B: At. Mol. Opt. Phys. 37, 763 (2004). 

[37] M. Ruggenthaler, D. Bauer, Phys. Rev. Lett. 102, 233001 (2009). 

[38] G. Tanner, K. Richter, and J.-M. Rost, Rev. Mod. Phys. 72, 497 (2000). 

[39] D. Hochstuhl, and M. Bonitz, proceedings of conferencne "PNGF IV", this issue (2009). 

[40] I. Kawata, H. Kono, A.D. Bandrauk, Phys. Rev. A 64 043411 (2001). 

[41] N. Suzuki, I. Kawata, and K. Yamashita, Chemical Physics 338, 348-353 (2007). 

[42] J.C. Light, I.P.Hamilton, and J.V. Lill, J. Chem. Phys. 82, 1400 (1985). 

[43] N.E. Dahlen, and R. van Leeuwen, J. Chem. Phys. 122, 164102 (2005). 

[44] N.E. Dahlen, A. Stan, and R. van Leeuwen, J. Phys: Conf. Ser. 35, 324 (2006); N.E. Dahlen, R. van 

Leeuwen, and A. Stan, J. Phys: Conf. Ser. 35, 340 (2006). 
[45] A. Stan, N.E. Dahlen, and R. van Leeuwen, J. Chem. Phys. 130, 224101 (2009). 
[46] J.C. Light, and T. Carrington Jr., Advances in Chemical Physics 114, 263 (2007) 

[47] The generalization to different numbers of basis functions per element is straightforward, and only slightly 

alters the matrix elements involved. 
[48] D.E. Manolopoulos and R.E. Wyatt, Chem. Phys. Lett. 152, 23 (1988). 

[49] The two-electron integrals u^jj are symmetric with respect to interchange of i <-> i <-> j and tupels 

[50] S. Bauch, and M. Bonitz, Phys. Rev. A 78, 043403 (2008). 
[51] W. Ku and A.G. Eguiluz, Phys. Rev. Lett. 89, 126401 (2002). 

[52] J.D. Alexander, C.R. Calvert, R.B. King, O. Kelly, L. Graham, W.A. Bryan, G.R.A.J. Nemeth, W.R. Newell, 
C.A. Froud, I.C.E. Turcu, E. Springate, I.D. Williams, and J.B. Greenwood, J. Phys. B: At. Mol. Opt. 
Phys. 42, 141004 (2009). 



